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AN  EFFICIENT  IMPEEMENTATION  OF  A 
BATCH-ORIENTED,  MULTITARGET,  MULTIDIMENSIONAL  ASSIGNMENT 
TRACKING  ALGORITHM  WITH  APPLICATION  TO  PASSIVE  SONAR 


I.  INTRODUCTION 


This  study  investigates  the  implementation  of  two  versions  of  a  batch-oriented, 
multidimensional  assigmnent  (MDA)  approach  for  tracking  contacts  with  long-duration 
crossings,  on  the  order  of  100  scans  for  a  single  dimensional  space.  The  environment 
being  simulated  is  that  of  the  passive  sonar  string  in  a  line  array  with  bearings-only 
information.  The  model  incorporates  multitarget  tracking  (MTT)  in  a  cluttered 
environment  with  missing  measurements.  Linear  programming  (LP)  relaxation  is  used  to 
solve  the  cost  assignment  matrix.  The  assignment  costs  are  calculated  via  the  Kalman 
filter  likelihood  function.  A  constraint  matrix  is  set  up  for  the  various  batch  lengths,  and 
the  assignment  problem  is  solved  via  an  LP  package.  This  study  concentrates  on  using 
efficient  techniques  to  eliminate  many  of  the  redundancies  in  generating  the  cost 
assignment  matrix,  and  a  single-scan,  forward-looking  filtering  method  is  developed  to 
further  reduce  the  number  of  generated  cost  evaluations. 

The  algorithm  is  generic  and  can  be  extended  to  the  A-dimcnsional  (ND) 
assignment  problem.  The  term  ND  refers  to  a  number  N  of  dimensions  D  where  the  first 
dimension  refers  to  the  tracker  state  estimates  and  covariance  from  the  previous  scan,  and 
( N  ~  1)  dimensions  refers  to  the  batch  length  in  time  for  the  subsequent  scans.  To 
minimize  computational  requirements,  a  suboptimal  version  of  the  algorithm  was 
developed  and  tested.  The  suboptimal  version  is  a  faster  implementation  of  the  ND 
assignment  algorithm,  and  results  have  been  processed  up  to  16D.  Simulated  results  have 
been  processed  up  to  4D  for  the  standard  implementation  of  the  algorithm.  The  efficient 
implementation  utilizes  (1)  a  two-dimensional  (2D)  LP  module  as  a  filtering  procedure  to 
extract  measurements  in  clutter  and  (2)  a  skip  factor  of  the  length  of  the  batch  in 
processing  the  ND  assignment.  The  standard  implementation  processes  the  ND 
assignment  and  keeps  only  the  first  scan’s  measurement  of  the  batch,  with  the  procedure 
repeated  for  consecutive  scans  via  a  sliding  batch. 
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Comparisons  are  performed  for  various  cases.  Test  data  are  generated  in  a 
cluttered  environment  via  simulation  and  processed  with  various  clutter  parameter  and 
batch  settings.  In  this  study,  only  simulations  with  linear  motion  model  targets  were 
considered,  and  the  test  data  simulate  detections  obtained  in  a  passive  sonar  environment 
for  target  crossings.  Monte  Carlo  runs  are  conducted  for  100  trials  for  the  various  cases. 
The  results  show  that  the  suboptimal  but  fast  implementation  tracked  correctly  through 
long-duration  crossings  in  a  much  larger  percentage  of  the  trials — over  90% — than  the 
standard  implementation — less  than  30%  of  the  trials. 

A.  BACKGROUND 

A  passive  line  array  in  a  typical  ocean  environment  detects  large  vessels  such  as 
tankers  and  freighters  at  a  long  distance,  so  that  contact  information  is  generally 
“bearings  only.”  When  two  of  these  distant  contacts  cross  in  bearing,  the  crossing 
duration  can  be  on  the  order  of  minutes.  During  the  crossing,  if  the  contacts  are  not  of 
approximately  equal  signal-to-noise  ratio  (SNR),  detections  will  be  assigned  to  the  higher 
SNR  contact,  and  the  lower  SNR  track  will  be  lost.  For  contacts  that  have  equal 
probability  of  detection  (Pd),  Willett  et  al.  [1]  have  recently  shown  that  the  estimated 
tracks  generated  with  a  standard  assignment  model  exhibit  a  repulsion  behavior  during 
the  crossing  that  violates  the  expected  behavior  of  these  large  contacts,  i.e.,  that  they  are 
known  to  travel  at  constant  course  and  speed  and  cannot  make  sudden  changes  because  of 
their  size.  This  repulsion  behavior  causes  problems  downstream  for  algorithms  that 
require  state  estimates  to  be  consistent  with  contact  trajectory. 

This  study  concentrates  on  utilizing  a  batch-style  algorithm  to  track  targets  with 
equal  Pd  in  a  long-duration  crossing.  MDA  techniques  (generally  batch-oriented)  are 
used  to  investigate  tracking  of  long-duration  crossings.  Willet  et  al.  [1]  explored  the  use 
of  batch-oriented  tracking  algorithms  for  the  linear  crossing  case.  The  algorithms 
discussed  in  this  thesis  can  be  generalized  to  N  dimensions.  Solving  this  A-dimensional 
batch  assignment  is  considered  to  be  “NP  hard”  (where  NP  means  nondeterministic 
polynomial  time)  for  the  case  in  which  the  dimensional  value  N  >  3;  therefore,  linear 
programming  relaxation  techniques  are  used.  An  off-the-shelf  software  package 
(LP  SOLVE)  [2],  [3]  is  used  to  solve  the  objective  function  for  the  assignment  problem. 


2 


Three  algorithms — a  modified  probabilistic  multi-hypothesis  tracker  (PMHT),  a  multi¬ 
hypothesis  tracker  (MHT),  and  an  MHT  with  rollout — were  investigated  and  tested  in  [1], 
and  those  simulations  and  test  results  are  the  foundation  for  the  simulated  results  in  this 
study. 

In  the  passive  sonar  environment,  the  contact  detections  from  a  line  array  are 
processed  in  sequence  by  a  beamformer  and  a  detector  and  subsequently  sent  to 
postprocessing  algorithms  such  as  tracking  and  classification  algorithms.  Tracks  are  then 
formed  for  the  postdetection  data  and  sent  to  other  postprocessing  algorithms  that  require 
tracks  as  inputs,  such  as  target  motion  analysis  (TMA)  algorithms  [4],  [5],  [6].  In  the 
scope  of  this  thesis,  simulations  are  developed  that  represent  the  postbeamformer  output 
of  a  detector  for  a  line  array  of  hydrophones. 

A  beamformer  [7],  [8]  is  essentially  a  transfonnation  from  the  hydrophone’s  or 
element’s  “time”  space  to  “bearing”  space  (i.e.,  beam  angle  or  conical)  space.  In  this 
thesis,  a  line  array  of  elements  is  considered.  A  beamformer  from  a  line  array  of 
elements  produces  a  fixed  number  of  beams,  based  on  directional  angles  (bearings),  by 
delaying  and  summing  the  time  series  from  each  element.  A  line  array  has  only  a  single 
dimension  for  its  observation  space,  which  is  its  conical  beam  space.  The  outputs  of  the 
beamformer  are  sent  to  a  detection  process  that  will  filter  and  nonnalize  the  data  to 
eliminate  the  unneeded  frequency  ranges  and  preserve  the  ranges  of  interest.  The 
detector  outputs  are  also  assumed  to  be  normalized  spatially  to  have  a  constant  noise 
background.  It  is  assumed  in  this  study  that  the  outputs  of  the  detector  are  peak-picked. 
A  peak-picked  value  is  any  measurement  in  the  detector  output  above  a  preselected 
threshold.  Because  the  beamformer  has  a  fixed  number  of  beams,  the  peak-picked  values 
from  the  detector  provide  only  discrete  measurements.  These  discrete  measurements  are 
interpolated  based  on  amplitude  or  other  criteria  to  provide  a  finer  estimate  of  the 
measurements  at  each  scan.  These  fine  measurements  will  be  sent  to  a  tracking 
algorithm. 

The  crossing  target  case  considered  is  the  long-duration  crossing  following  a 
linear  trajectory.  Test  results  are  processed  via  Monte  Carlo  simulation  for  various 
clutter  levels  and  various  batch  lengths. 
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B. 


REVIEW  OF  BATCH-ORIENTED  MTT  ALGORITHMS 


Recent  developments  in  MTT  and  data  fusion  technologies  are  all  pointing  to  the 
need  to  include  batch  frames  of  data  in  order  to  achieve  optimal  performance.  The 
current  batch-oriented  MTT  algorithms  include  four  classes — the  PMHT,  the  MHT,  the 
batch-oriented  joint  probabilistic  data  association  (JPDA)  tracker,  and  the  MDA  tracker. 
The  PMHT  developed  by  Streit  and  Luginbuhl  [9],  [10]  is  one  batch-oriented  algorithm 
that  has  shown  promise  for  passive  sonar  and  other  applications.  However,  the  PMHT 
does  not  directly  produce  accurate  error  covariance  matrices  for  the  track  estimates;  a 
separate  computation  is  required  [1].  The  MHT,  as  originally  proposed  by  Reid  [11],  has 
a  batch-oriented  framework  where  tracks  are  constructed  based  on  enumeration  of  all 
possible  measurement-to-track  association  hypotheses  along  a  batch.  As  time  evolves, 
the  number  of  tracks  grows  exponentially,  based  on  new  measurement  arrivals  with  each 
scan.  For  practical  implementation,  current  MHT  implementations  are  suboptimal  [12] 
because  of  computational  processing  and  memory  limitations  associated  with  using  fixed- 
size  batch  lengths.  To  limit  the  exponential  growth  of  the  number  of  tracks,  even  with  a 
fixed  batch  length,  ad  hoc  logic  is  used  to  prune  tracks  that  are  infeasible  and  meet 
certain  criteria  at  a  subsequent  scan.  This  type  of  deferred  decision  logic  can  only  be 
applied  when  using  a  batch-style  tracker.  The  JPDA  algorithm  [13,  pp.  310-319],  [14]  is 
a  true  multitarget  tracking  algorithm  that  produces  consistent  state  estimates  (i.e., 
accurate  state  and  covariance  estimates),  but  batch  extensions  [15],  [16],  [17]  are  still 
very  limited  (batch  lengths  of  less  than  three  scans)  in  the  current  state  of  the  art. 
Further,  the  current  version  of  the  algorithm  is  prone  to  track  segmentation,  because  of 
the  track  coalescence  effect  of  PDA  style  trackers  [18].  In  PDA-style  trackers,  the  tracks 
tend  to  merge,  i.e.,  coalesce,  as  they  cross,  which  requires  initialization  of  a  new  track  to 
replace  the  segment  not  tracked  past  the  crossing.  The  MDA  algorithm  [19],  [20]  utilizes 
an  optimization  framework,  and  an  enumeration  of  all  of  the  possible  measurement-to- 
track  costs  is  calculated  along  the  batch.  Cost  minimization  is  perfonned  via  several 
techniques.  Lagrangian  relaxation  is  used  in  [21],  whereas  the  interior  point  linear 
programming  is  used  in  [22],  In  [2],  a  mixed  integer  linear  programming  (LP)  relaxation 
method  is  utilized  to  minimize  the  assignment  cost  matrix  with  an  open  source  solver, 
LP  SOLVE  [3].  The  MDA  framework  is  based  on  0-1  integer  assignments  as  originally 
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proposed  in  [23]  via  0-1  integer  linear  programming,  where  tracks  are  given  discrete 
assignments  with  one  measurement  per  scan  and  are  not  allowed  to  share  measurements, 
based  on  explicit  constraints  placed  on  the  optimization.  An  LP-based  method  to  produce 
mixed  integer  assignments  is  explored  in  [24].  The  current  work  will  investigate  an 
implementation  of  the  MDA  algorithm  that  uses  linear  programming  and  incorporates 
longer  batch  lengths  into  this  framework.  The  motivation  for  longer  batch  lengths  is  the 
premise  that  extending  the  batch  length  to  cover  the  tracks  for  periods  before  and  after 
the  crossing  period  will  aid  in  the  long-duration  crossing  target  problem  by  providing 
enough  data  on  the  contact  positions  to  better  predict  the  correct  trajectories.  This 
algorithm  will  be  analyzed  with  simulated  results  related  to  the  passive  sonar,  long- 
duration  crossing  tracks  problem.  A  suboptimal  method  is  also  developed  to  reduce  the 
computational  load  of  the  algorithm — by  reducing  the  number  of  cost  evaluations — to 
further  increase  the  batch  length. 

This  thesis  is  organized  as  follows.  In  Chapter  II,  the  two  algorithms  are 
described  and  the  motion  model  and  measurement  model  assumptions  are  fonnulated. 
An  implementation  of  the  MDA  algorithm  is  presented  as  Algorithm  I.  The  2D-ND 
MDA  algorithm,  which  is  a  more  efficient — and  therefore  faster — version,  is  presented  as 
Algorithm  II.  In  Chapter  III,  which  presents  the  results,  the  two  algorithms  are  applied  to 
the  long-duration  crossing  target  problem  for  a  passive  sonar  line  array.  Simulation 
results  are  presented  for  various  test  cases. 
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II.  ALGORITHM  DESCRIPTION 


A.  MOTION  MODEL  ASSUMPTIONS 


The  true  positional  state  vector  for  a  contact  at  time  scan  4+1  is  given  by 


x(Oi)  =  A 


40 

40 


(i) 


where  t  is  the  time  of  the  scan,  s  is  the  scan  index  of  the  original  data,  and  A  is  the  state 
transition  matrix  of  the  system. 


A  = 


1 

0 


(2) 


where  T  is  the  scan  period. 


B.  MEASUREMENT  MODEL  ASSUMPTIONS 


The  measurement  model  consists  of  received  measurements  normally  distributed 
about  the  fractional  beam  space  x. 

40  =40+40°;,  (3) 

R  =  °r\  (4) 

where  z(4)  is  the  linear  measurement  position  at  time  scan  4,  w  is  the  white  Gaussian 
random  noise  at  time  scan  4,  ar  is  the  standard  deviation  of  the  measurement  noise,  and  R 
is  the  variance  of  the  measurement  noise. 

The  clutter  measurements  are  assumed  to  have  a  Poisson  random  variable 
distribution  with  spatial  clutter  density  parameter  X.  The  contact’s  probability  of 
detection  PD  is  fixed,  and  spatial  clutter  density  parameter  X  is  varied  for  the  simulated 
results. 
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C.  ALGORITHM  I:  MDA  ALGORITHM 

1.  Overview  of  the  MDA  Algorithm 

An  overview  of  the  MDA  algorithm  is  provided  as  a  flowchart  in  Figure  1.  The 
flow  of  the  modules  is  based  on  [19]  and  [20],  except  for  the  LP  Assign  section,  which  is 
based  here  on  LP  relaxation  instead  of  the  Lagrangian  relaxation  technique.  In  this  study, 
it  is  assumed  that  the  tracks  are  already  initialized,  where  prior  track  estimates  are  known. 
Based  on  the  simulated  data,  the  tracks  are  initiated  on  the  first  scan  and  the  algorithm 
processes  the  subsequent  scans  based  on  the  initiated  tracks  via  a  sliding  batch  length  ( N 
-  1).  No  track  management  functions,  such  as  dropping,  pruning,  or  merging  tracks,  are 
considered  in  this  study.  No  new  tracks  will  be  initiated  subsequently  after  the  first  scan. 
Much  of  the  processing  for  an  MDA-style  algorithm  takes  place  in  the  GenCost  and 
LP  Assign  sections  of  Figure  1.  The  GenCost  section  computes  the  Kalman  filter 
likelihoods  for  every  possible  combination  of  measurement-to-track  associations  per  scan 
for  the  entire  (A-  1)  batch  length.  It  also  accounts  for  missing  measurements  per  scan 
for  the  entire  batch  length  per  track.  This  study  concentrates  largely  on  the  GenCost 
section  of  the  algorithm  and  eliminates  many  of  the  deficiencies  of  this  module.  The 
LP  Assign  section  performs  the  actual  measurement-to-track  assignment  via  LP 
relaxation.  This  assignment  can  be  performed  via  various  methods  such  as  Lagrangian 
relaxation,  interior  point  linear  programming,  and  integer  linear  programming  [2],  [21]- 
[24].  This  study  employs  linear  programming  relaxation  via  an  available  open  source 
solver,  LP  SOLVE  [3]. 
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Figure  1.  Overview  of  MDA  tracking  algorithm. 


2.  Track  Initialization 

The  track  initialization  function  initializes  a  track  based  on  known  input 
parameters  such  as  initial  target  (i.e.,  contact)  position  )  and  rate  x(t  )  .  As  part  of 

this  study,  the  tracks  are  initialized  on  the  initial  simulated  state  estimates  based  on  the 
motion  model  assumptions.  The  rate  term  x(tSo)  is  set  to  0  for  the  initialized  track’s 

state.  This  function  is  perfonned  at  the  first  scan. 
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3. 


Cost  Generation 


The  generation  of  the  cost  evaluations  used  for  the  objective  function  of  the  MDA 
algorithm  was  analyzed  for  efficiency  in  this  research,  and  the  resulting  GenCost  module 
is  described  here  in  detail.  Cost  evaluations  are  calculated  based  on  all  possible 
measurement-to-track  assignment  hypotheses  including  the  missing  measurement  along  a 
batch.  The  cost  evaluations  are  obtained  by  the  Kalman  filter  negative  log  likelihood 
calculation.  For  each  track,  a  missing  measurement  is  also  accounted  for  in  every  scan  of 
the  batch.  Note  that  the  following  prediction  equations  are  used  only  to  compute  the 
likelihoods.  At  time  tn,  the  start  of  each  batch,  the  state  estimates  and  covariance  for  each 
track  from  the  prior  scan  are  used  and  the  negative  log  likelihoods  are  cumulated  along 
the  batch.  These  prior  state  estimates  and  covariances  are  from  the  Kalman  filter  update 
section  of  this  thesis. 

Equations  (5)  through  (13)  are  from  standard  Kalman  filter  theory.  The  Kalman 
filter  time  update  comprises  the  state  prediction,  given  by  Equation  (5)  and  the  state 
prediction  covariance,  given  by  Equation  (6). 


x(t„+1)  =  Ax(0, 

(5) 

P(f„+1)  =  AP(OA'+Q, 

(6) 

where  the  index  n  is  the  scan  index  in  the  batch,  P  is  the  predicted  state  prediction 
covariance,  P  is  the  state  prediction  covariance  in  current  time,  and  Q  is  the  covariance  of 
the  discrete -time  process  as  defined  in  Equation  (7). 

P(0)  is  the  initial  covariance  matrix  set  at 

P(0)  = 

where  the  Kalman  filter’s  initial  function  value  variance  a  f2  is  a  parameter  and  is  the 
same  for  all  tracks. 
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The  covariance  matrix  Q  is  based  on  a  discrete  white  noise  acceleration  model 
from  [25,  p.  274]  and  given  by  Equation  (7): 


Q 


rj-i  4 

T 

T3 

T 


T3 

T 

T2 


(V) 


The  measurement  prediction  is  given  by 

z(tn+l)  =  Cx(tn+l),  (8) 

where  C  =  [1  0]  is  the  measurement  matrix,  and  aq  is  the  standard  deviation  used  in  the 
process  noise  model. 

The  Kalman  filter  measurement  update  is  computed  using  the  innovation 
covariance  estimate  S(fi+I ),  filter  gain  estimate  K(t„+i),  measurement  residual  v{tn+l), 
updated  state  estimate  \(/i+] ),  and  the  updated  state  covariance  P(/,+l ) ,  given  by 


Equations  (9)  through  (13),  respectively. 

S(tn+l)  =  CP(tll+l)C'  +  R,  (9) 

K(t/1+1)  =  P(t„+1)C'S(t/I+ir1,  (10) 

v(tn+l)  =  z(tn+l)~z(tn+l),  (11) 

x(t,.+1)  =  x(t„+1)  +  K (tn+1)v(tnJ ,  (12) 

ntM)  =  P(^)-K(t„+1)S(t„+I)K(t„+1)' .  (13) 


The  cost  calculation  [24],  [26]  is  perfonned  via  the  following  function  for  the 
Kalman  filter  likelihood  A: 


A«„,  >  =  A(r, )  +  )-' V«„, )  +  In 


jfel  ^^m)' 


(14) 


To  process  a  missing  measurement,  a  penalty  is  imposed  where 
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A('„+I)  =  A(0-ln(l-^)- 


(15) 


Note  that  Ae  is  the  assumed  clutter  density  used  by  the  tracking  algorithm  in 
Equation  (14).  The  clutter  density  parameter  differs  from  A  in  the  measurement  model 
section  only  in  that  Ae  is  a  fixed  parameter  whereas  A  varies  with  the  simulation 

scenario.  Note  that,  in  the  results  section,  the  case  is  examined  where  Ae  equals  A,  i.e., 

the  matched  case.  If  a  missing  measurement  is  being  processed,  only  the  Kalman  time 
update  Equations  (5)  through  (8)  and  (15)  are  performed.  This  set  of  equations — 
Equations  (5)  through  (15) — is  used  to  calculate  the  Kalman  filter  likelihoods  A  that 
serve  as  the  costs  forming  the  objective  function  for  the  MDA  algorithm.  The  basic 
computational  unit  for  an  MDA  algorithm  is  a  single  Kalman  filter  update  used  to 
calculate  the  negative  log  likelihoods.  The  details  of  the  negative  log  likelihood  function 
are  described  in  [24]  and  [26].  Reference  [21]  provides  further  details  on  the  cost 
calculation  methodology  used  for  the  MDA  algorithm.  The  calculations  of  Equations  (5) 
and  (6)  and  (8)  through  (14)  are  performed  for  every  track  at  each  scan  along  the  batch 
for  all  possible  combinations  of  measurement-to-track  associations.  For  missing 
measurements,  Equations  (5),  (6),  (8),  and  (15)  are  calculated.  Note  that  Equations  (14) 
and  (15)  are  used  for  processing  actual  measurements  or  missing  measurements 
respectively.  A  single  likelihood  represents  a  possible  trajectory  for  a  track  along  a 
batch.  Enumeration  of  all  measurement-to-track  (including  missing  measurement) 
likelihoods  needs  to  be  calculated  per  track  for  the  entire  batch  length.  The  number  of 
costs  calculated  per  track  is  given  by  Equation  (16)  for  a  fixed  number  of  measurements 
per  scan  along  a  batch. 

N- 1 

Number  of  costs  (c)  =  MN~X  or  ,  (16) 

n= 1 

where  M  is  the  fixed  number  of  measurements  per  scan  n,  including  the  missing 
measurement,  and  ( N~  1)  is  the  batch  length.  The  total  number  of  costs  is  given  by 

N- 1 

Total  number  of  costs  =  Mn~'L  or  E  ,  (17) 

n= 1 

where  L  is  the  number  of  tracks. 
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The  total  number  of  Kalman  filters  (KF)  to  process  for  all  tracks  in  a  batch  is 
given  by 

AM 

Total  number  of  KF  =  Mn~1L(N-  1)  or  L-(N-\)-\[Mn  .  (18) 

n= 1 

The  Kalman  filters  used  to  compute  the  likelihoods  are  the  basic  unit  as  far  as 
computational  costs  for  a  MDA-style  algorithm,  since  it  must  be  performed  numerous 
times  to  evaluate  all  possible  combinations  of  measurement-to-track  costs.  Note  that 
Equations  (16)  through  (18)  apply  only  to  cases  with  a  fixed  number  of  measurements  per 
scan  and  can  be  used  for  rough  algorithm  loading  calculations. 

If  the  number  of  measurements  changes  from  scan  to  scan,  as  normally  happens, 
the  cost  tree  and  the  number  of  Kalman  filters  are  calculated  using  generic  equations. 
The  number  of  costs  for  any  number  of  measurements  and  the  number  of  Kalman  filters 
per  track  are  calculated  using  Equations  (19)  and  (20),  respectively. 

AM 

Number  of  costs  (c)  =  II'”..  <19) 

n= 1 

N-l 

Number  of  KF  updates  =  (iV-l)]^[/n„  ,  (20) 

n= 1 

where  m  is  the  varying  number  of  measurements  for  scan  n  from  1  to  (N~  1). 

Consider  an  example  of  two  tracks  with  a  batch  length  of  (N  -  1)  of  three  scans 
with  three  measurements  in  each  scan,  as  shown  in  Figure  2.  This  is  a  4D  Assignment 
with  a  batch  length  of  three.  The  measurements  include  the  missing  measurement,  where 
there  are  only  two  actual  measurements.  The  topology  on  the  left  represents  the  track 
paths  for  Track  1  and  on  the  right  for  Track  2.  The  measurements  are  labeled  from  one  to 
nine  for  the  three  measurements  and  three  scans  of  data  as  shown.  Note  that 
measurements  1,  4,  and  7  are  labeled  as  missing  measurements.  From  Equation  (16), 
there  are  3'  or  27  possible  paths  for  this  specific  example  per  track.  As  mentioned  in 
Equation  (18),  the  total  number  of  Kalman  filter  updates  to  process  is  Mn~1L(N-  1)  .  In 
this  case,  it  would  require  162  Kalman  filters  to  generate  a  total  of  54  cost  values  c  in  a 
single  batch  for  the  two  tracks. 
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It  is  important  to  note  that  there  is  no  gating  procedure  in  Algorithm  I.  In  the 
MDA  procedure,  as  outlined  in  [19]  and  [21],  there  is  no  gating  involved,  with  every 
measurement-to-track  hypotheses  calculated,  as  is  the  case  for  Algorithm  I.  To  improve 
processing  speed  of  MDA-style  algorithms  gating  can  be  perfonned  as  shown  in  [27], 
using  a  clustering  technique  along  the  batch.  Any  type  of  gating  used  makes  the 
association  problem  suboptimal,  especially  when  processing  increasing  time  depths.  A 
suboptimal  gating  technique  based  on  a  2D-ND  MDA  algorithm  is  proposed  for  faster 
processing.  The  gating  procedure  is  described  in  Section  D.3  for  Algorithm  II. 


Scan  3 

Scan  2 

Scan  1 

T  racks 

Figure  2.  Example  of  4D  assignment  case  cost  mapping  with  two  tracks  and  three 
measurements,  including  the  missing  measurement. 


In  Figure  3,  the  unraveled  track  paths  are  shown  for  Track  1;  the  measurement 
indices  in  this  figure  are  labeled  as  in  Figure  2.  Notice  the  tree  structure  of  the  different 
paths  from  the  track  to  the  end  of  the  batch.  The  main  trunk  of  the  tree  is  at  the  bottom 
where  the  tracks  start  from  the  first  scan,  and  the  canopy,  fringe  branches  are  towards  the 
last  scan. 
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Figure  3.  Example  of  4D  cost  tree  mapping  for  Track  1 . 


Table  1  provides  an  enumeration  of  the  binary  value  indexes  p  for  obtaining  the 
cost  evaluation  tree  depicted  in  Figure  3.  As  previously  mentioned,  for  the  4D  example, 
27  cost  terms  c  are  calculated  per  track.  The  top  row  represents  the  numbered 
measurements  depicted  in  Figure  2.  The  leftmost  column  represents  the  associated 
binary  decision  values  p  to  the  costs  c.  The  four  subscripted  numerals  to  the  decision 
term  p  represent  an  N-  tuple  path  along  the  batch  for  a  given  track.  In  Table  1,  the  first 
tenn  in  the  subscript  represents  Track  1  ’s  costs  and  associated  decision  terms  p,  or  the 
index  of  the  first  dimension.  The  second  term  represents  the  measurement  index  in  the 
first  scan,  the  third  term  represents  the  measurement  index  in  the  second  scan,  and  so  on. 
A  Kalman  filter  calculation  is  represented  by  a  value  of  1  in  Table  1,  and  for  this  4D 
example  there  are  8 1  filters  used  per  track. 

These  costs  enforce  the  underlying  assignment  model  where  the  constraints 
imposed  require  that  only  one  measurement  is  associated  with  a  track  per  scan  along  a 
batch  and  that  multiple  tracks  can  be  associated  with  the  missing  measurement  at  the 
same  scan.  This  optimization  procedure  is  described  for  Algorithm  I  in  Section  C.5. 


15 


In  Table  2,  the  number  of  calculated  costs  c  for  a  single  track  for  various  N 
dimensional  assignments  from  2D  to  10D  assignment  and  number  of  measurements  from 
1  to  10  per  scan  are  provided. 


\I\M 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

3 

1 

4 

9 

16 

25 

36 

49 

64 

81 

100 

4 

1 

8 

27 

64 

125 

216 

343 

512 

729 

1000 

5 

1 

16 

81 

256 

625 

1296 

2401 

4096 

6561 

10000 

6 

1 

32 

243 

1024 

3125 

7776 

16807 

32768 

59049 

100000 

7 

1 

64 

729 

4096 

15625 

46656 

117649 

262144 

531441 

1000000 

8 

1 

128 

2187 

16384 

78125 

279936 

823543 

2097152 

4782969 10000000 

9 

1 

256 

6561 

65536 

390625 

1679616 

5764801 1677721643046721 

1E+08 

10 

1 

512 

19683 

262144  1953125  10077696  40353607  1.34E+08  3.87E+08 

1E+09 

Table  2. 

Number  of  /V-tuple  cost  evaluations  c 

for  various  dimensions  N  with  a  fixed 

number  of  measurements  Mper  scan. 


Table  3  provides  the  number  of  Kalman  filters  needed  to  generate  the  costs  in 
Table  2. 


N\M 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

3 

2 

8 

18 

32 

50 

72 

98 

128 

162 

200 

4 

3 

24 

81 

192 

375 

648 

1029 

1536 

2187 

3000 

5 

4 

64 

324 

1024 

2500 

5184 

9604 

16384 

26244 

40000 

6 

5 

160 

1215 

5120 

15625 

38880 

84035 

163840 

295245 

500000 

7 

6 

384 

4374 

24576 

93750 

279936 

705894 

1572864 

3188646 

6000000 

8 

7 

896 

15309 

114688 

546875 

1959552 

5764801  14680064  33480783  70000000 

9 

8 

2048 

52488 

524288 

3125000  13436928  46118408  1.34E+08  3.44E+08 

8E+08 

10 

9 

4608  177147  2359296  17578125  90699264  3.63E+08  1.21E+09  3.49E+09 

9E+09 

Table  3.  Number  of  Kalman  filters  needed  to  generate  cost  c  for  various 

dimensions  N  with  a  fixed  number  of  measurements  Mper  scan. 
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As  noted  in  [27]  and  [28],  the  costs  calculated  using  the  GenCost  function,  or 
generating  the  Kalman  filter  negative  log  likelihood  A,  take  up  about  95%  of  the 
computational  requirement  of  an  MDA-style  of  algorithm.  In  Table  3  for  example,  with 
nine  measurements  and  7D  assignment,  over  three  million  Kalman  filter  updates  need  to 
be  processed  for  a  single  track  to  generate  all  the  costs  in  the  batch.  As  the  number  of 
tracks  increase,  the  cost  calculations  and  the  number  of  Kalman  filter  updates  to  process 
increase  linearly  with  the  numbers  from  Tables  2  and  3,  respectively.  If  the  cost 
generation  function  is  taking  up  much  of  the  computational  resources,  then  it  is  beneficial 
to  reduce  the  number  of  overall  Kalman  filter  calculations  along  a  batch  per  track.  Some 
of  the  cost  generation  inefficiencies  can  be  eliminated  because  the  Kalman  filter  negative 
log  likelihoods  have  previously  been  calculated  for  the  initial  scans  of  the  batch.  Storing 
a  select  set  of  previously  processed  cumulative  likelihoods,  prior  track  filter  state 
estimates,  and  covariances  makes  these  values  available  for  reuse  later  in  the  processing, 
which  results  in  significant  computational  savings  without  any  effect  on  the  overall 
algorithm.  The  circled  sections  in  Table  4  show  the  costs  that  are  part  of  the  main  trunk 
of  the  tree  that  can  be  reused  to  calculate  the  canopy  branches  of  the  cost  tree  from  Figure 
3.  It  is  observed  that  the  number  of  links  in  the  cost  tree  corresponds  to  the  number  of 
Kalman  filters  to  process.  The  main  trunk  sections  are  early  in  the  batch,  and  the  fringe 
branches,  representing  the  number  of  cost  evaluations,  are  toward  the  end  of  the  batch. 
By  using  this  technique,  the  number  of  Kalman  filter  updates  to  calculate  in  the  4D 
example  decreases  from  162  to  78  filters  for  the  two  example  tracks,  which  is  a 
significant  savings  without  any  impact  on  the  integrity  of  the  algorithm.  Note  that  the 
number  of  costs  calculated  for  the  objective  function  is  still  the  same.  The  reduction  in 
the  number  of  Kalman  filters  provides  the  computational  savings  to  process  the  costs. 

The  equation  used  to  calculate  the  number  of  reduced  Kalman  filter  (KF)  updates 
for  a  fixed  number  of  measurements  M  in  a  batch  is  given  by  the  following  special  case 
of  the  geometric  series  and  generation  function,  where  M  >  1 : 


JV-l 


Number  of  reduced  KF  =  y'  M " 


n= 1 


mn -m 

M- 1 


(21) 
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Even  when  processing  a  missing  measurement,  a  Kalman  filter  time  update  must 
be  performed  to  calculate  the  state  estimates  for  the  next  scan.  The  number  of  Kalman 
time  updates  to  process  missing  measurements  is  given  by 


JV— 1 

Number  of  reduced  KF  time  updates  =  2>""' 

n= 1 


-1 

M- 1 


(22) 


The  total  number  of  full  Kalman  filters  to  process  real  measurements  along  a 
batch  is  given  by 


JV— 1 

Number  of  reduced  full  KF  =  y = 

n= 1 


mn  -Mw 

M  -1 


-1  = 


-1.  (23) 


The  generic  equation  to  calculate  the  number  of  reduced  Kalman  filters  for  any 
number  of  measurements  in  a  batch,  for  a  single  track  is  given  by 

N  k-l 

Number  of  reduced  KF  =  E  IT”.  •  <24) 

k=2  n= 1 


where  mn  is  the  number  of  measurements  for  scan  n  from  I  to  (k  -  1 ),  and  k  is  the 
dimensional  assignment  from  2  to  N. 

Table  5  shows  the  number  of  Kalman  filters  required  using  the  efficient 
methodology  used  in  this  study  for  N  dimensions  from  2D  to  10D  and  a  fixed  number  of 
measurements  per  scan,  from  1  to  10.  From  Table  3,  as  was  mentioned  previously  in  an 
example,  with  nine  measurements  and  7D  assignment,  over  three  million  Kalman  filters 
needed  to  be  processed  for  a  single  track.  In  Table  5,  with  the  reduced  number  of 
Kalman  filters  with  nine  measurements  and  7D  assignment,  the  number  of  Kalman  filters 
dramatically  decreases  to  less  than  600,000 — a  more  than  fivefold  savings. 

Table  6  shows  the  savings  in  computation  achieved  with  the  reduced  set  of 
Kalman  filters  based  on  the  efficient  storage  methodology.  This  table  gives  the  ratio  of 
the  number  of  Kalman  filters  for  each  case  in  Table  3 — the  full  set — to  the  number  of 
Kalman  filters  required  for  the  same  case  in  Table  5 — the  reduced  set.  As  is  expected,  as 
the  N  dimensional  size  and  the  number  of  measurements  per  scan  increase,  the 
computational  savings  significantly  increase.  In  the  10D  case  with  10  measurements  per 
scan,  there  is  an  eightfold  decrease  in  the  number  of  Kalman  filters  computed. 
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N\M 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

3 

2 

6 

12 

20 

30 

42 

56 

72 

90 

110 

4 

3 

14 

39 

84 

155 

258 

399 

584 

819 

1110 

5 

4 

30 

120 

340 

780 

1554 

2800 

4680 

7380 

11110 

6 

5 

62 

363 

1364 

3905 

9330 

19607 

37448 

66429 

111110 

7 

6 

126 

1092 

5460 

19530 

55986 

137256 

299592 

597870 

1111110 

8 

7 

254 

3279 

21844 

97655 

335922 

960799 

2396744 

5380839 

11111110 

9 

8 

510 

9840 

87380 

488280 

2015538 

6725600  19173960  48427560 

1.11E+08 

10 

9 

1022 

29523 

349524  2441405  12093234  47079207  1.53E+08  4.36E+08 

1.11E+09 

Table  5.  Number  of  Kalman  filters  needed  in  reduced  set  where  previously 

calculated  costs  c  and  filter  outputs  are  stored  for  various  dimensions  N 
with  a  fixed  number  of  measurements  M  per  scan. 


N\M 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

3 

1 

1.33 

1.50 

1.60 

1.67 

1.71 

1.75 

1.78 

1.80 

1.82 

4 

1 

1.71 

2.08 

2.29 

2.42 

2.51 

2.58 

2.63 

2.67 

2.70 

5 

1 

2.13 

2.70 

3.01 

3.21 

3.34 

3.43 

3.50 

3.56 

3.60 

6 

1 

2.58 

3.35 

3.75 

4.00 

4.17 

4.29 

4.38 

4.44 

4.50 

7 

1 

3.05 

4.01 

4.50 

4.80 

5.00 

5.14 

5.25 

5.33 

5.40 

8 

1 

3.53 

4.67 

5.25 

5.60 

5.83 

6.00 

6.13 

6.22 

6.30 

9 

1 

4.02 

5.33 

6.00 

6.40 

6.67 

6.86 

7.00 

7.11 

7.20 

10 

1 

4.51 

6.00 

6.75 

7.20 

7.50 

7.71 

7.88 

8.00 

8.10 

Table  6.  Ratio  of  number  of  Kalman  filters  calculated  for  full  set  (per  Table  3)  vs 
reduced  set  (per  Table  5).  In  both  sets,  N  is  the  dimensional  size,  and  there 
is  a  fixed  number  of  measurements  Mper  scan. 


The  storage  of  previously  calculated  costs  c  and  Kalman  filter  outputs  is  used  in 
this  version  of  MDA  Algorithm  I  and  in  the  suboptimal  but  faster  Algorithm  II. 
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4. 


Linear  Programming  Overview 


Once  the  costs  c  are  obtained,  based  on  an  enumeration  of  A-tuple  measurement- 
to-track  paths  along  a  batch,  the  assignment  of  these  costs-to-track  can  be  performed  via 
a  variety  of  methods.  In  [21],  an  efficient  near-optimal  algorithm  is  developed  using  a 
relaxation  algorithm  based  on  Lagrangian  multipliers,  with  a  modified  Auction  algorithm 
[29].  Other  efficient  methods  are  also  currently  available  for  solving  the  objective 
function  for  global  constraint  optimization  problems  [2].  Some  of  the  other  solutions 
were  mentioned  in  Chapter  I  in  the  review  of  batch-oriented  MTT  algorithms.  There  is 
also  a  variety  of  global  optimization  solvers  commercially  available  or  as  open  source 
software  packages.  As  part  of  this  study,  an  open  source  package  called  LPSOLVE  [2], 
[3]  is  used  to  perform  the  linear  global  constraint  optimization.  LP  SOLVE  utilizes  a 
mixed  integer  linear  programming  relaxation  technique  via  the  primal-dual  method. 
Linear  programming  (LP)  was  originally  developed  by  Dantzig  [30]  in  1947  and  was 
initially  used  by  the  U.S.  Air  Force.  Currently,  LP  is  used  in  a  variety  of  fields  to  solve 
the  linear  global  constraint  optimization  problem.  LP  sets  up  the  optimization 
formulation  in  the  fonn  of  Ax  =  B,  where  A  contains  the  associated  V- tuple  costs  c,  x 
contains  the  binary  decision  variable  p,  and  B  contains  the  associated  constraint  terms  as 
outlined  in  this  study.  A  detailed  description  of  LP  is  beyond  the  scope  of  this  thesis,  but 
the  reader  is  referred  to  several  texts  [30],  [31],  [32]  on  this  well-researched  and  much- 
applied  method.  The  ND  assignment  operation  is  performed  in  the  LP  Assign  module  of 
Figure  1,  and  objective  formulation  is  outlined  in  the  next  section. 

5.  ND  Assignment 

The  ND  assignment  refers  to  a  batch  length  of  (A  -  1)  scans  where  the  first 
dimension  represents  the  prior  state  estimates  and  covariance  of  a  given  track.  For 
example,  a  16D  assignment,  or  a  batch  length  of  15  scans,  is  processed  for  each  target. 
The  ND  assignment  formulation  follows  Equations  (25)  and  (26). 
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The  objective  function  is  given  by 


min 

Pin- "in 


ml  m2  mn 

ZZ-Z 

i\=u\  h  =u2  in=un 


c 

iVT-in 


(25) 


where  c  is  the  /V-tuple  Kalman  filter  negative  log  likelihood,  subject  to  the  set  of 
constraints  defined  as 


m2  m2  mn 

ZZ-Z/v 

i2=u2  z3  =«3  in  =un 


=  1  for  i  =  1, 2,  •  •  •  ml 


m]  m3 

ZZ-Zp 

zj  =«!  Z3=tt3  in  =un 


for  i  =1,2 ,•••«, 

2  z 


(26) 


A7t]  aw2  m;(  —1 

Z  £•••  Z  =1  for  z'„  =^--mn, 

i\=U\  =u2  in-\=Un 

where  /«„  is  the  number  of  measurements  for  scan  n  from  1  to  (N—  1).  The  measurement 
index  i  ranges  from  un  to  mn,  where  u„  =  0,  and  n  is  the  scan  index.  Note  that  (un  =  0)  is 
the  missing  measurement.  The  values  of  c‘i|/2...,  are  the  /V-tuple  costs.  These  costs  are 

the  cumulative  Kalman  filter  negative  log  likelihoods  for  the  batch,  A  ,  as  obtained  in  the 
GenCost  function.  The  costs  are  minimized  in  Equation  (25)  subject  to  the  constraints  in 
Equation  (26).  The  constraints  impose  the  requirement  that  at  most  one  measurement  is 
associated  with  a  track  per  scan  in  a  batch  and  multiple  tracks  can  be  associated  with  the 
missing  measurement  at  the  same  scan.  In  other  words,  each  measurement-to-track 
assignment  is  exclusive  per  scan  in  the  batch.  It  is  noted  in  [27]  that  the  actual 
measurement-to-track  assignment  is  a  very  efficient  process  and  only  takes  5%  of  the 
computational  resources  for  the  algorithm.  Once  the  measurement-to-track  assignments 
are  obtained,  the  next  step,  as  shown  in  Figure  1,  is  to  update  the  track’s  state  estimate  via 
a  standard  Kalman  filter  as  described  in  the  next  section. 
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6. 


Kalman  Filter  Update 


Once  the  measurement-to-track  assignment  is  performed,  a  Kalman  filter  is  used 
to  update  the  state  matrix  for  the  next  scan  ts.  The  filter  process  is  performed  via  the 
Kalman  filter  state  space  Equations  (5)  through  (13),  calculated  for  each  track  at  current 
scan  ts  with  measurement  z(t\)  obtained  for  the  track  via  the  (LP_Assign)  assignment 
function  for  the  first  scan  of  the  batch.  Note  that,  in  Equations  (5)  through  (13),  the  scan 
index  n  represents  the  scan  index  in  a  batch,  which  is  replaced  in  this  step  with  s,  the  scan 
epoch  index  for  the  total  number  of  scans. 

Note  that  an  assignment  of  a  missing  measurement  only  performs  the  Kalman 
time  update,  Equations  (5)  through  (8).  The  filter  is  processed,  and  the  cycle  repeats  for 
the  next  scan  of  data  according  to  Figure  1. 

D.  ALGORITHM  II:  A  FASTER,  SUBOPTIMAL  VERSION  VIA  2D-ND  MDA 

1,  Overview  of  Algorithm  II 

The  primary  goal  of  this  research  is  to  extend  the  batch  length  size  of  the  tracker 
with  the  premise  that  it  will  aid  in  the  long-duration  target  crossings  problem.  In 
Algorithm  I,  it  is  shown  that  by  increasing  the  dimensional  size,  the  number  of  Kalman 
filters  to  process  the  cost  evaluations  increases  exponentially.  In  [27],  it  is  noted  that 
even  though  the  assignment  function  of  the  MDA  is  the  most  important  part  of  the 
algorithm,  it  accounts  for  only  5%  of  the  computational  resources.  About  95%  of  the 
computational  load  is  due  to  the  generation  of  the  cost  evaluations,  to  form  the  objective 
function,  given  by  Equation  (25).  In  Algorithm  I,  the  number  of  cost  evaluations 
remained  the  optimal  size  with  a  full  enumeration  of  all  A-tuple  pairwise  measurement- 
to-track  decisions  along  the  batch  were  calculated.  In  [27],  a  fast  MDA-style  algorithm  is 
developed  by  using  a  clustering  process  to  reduce  the  number  of  cost  evaluations.  This  is 
a  type  of  gating  procedure  along  a  batch  based  on  a  clustering  approach.  As  part  of  this 
research,  in  Algorithm  II,  a  method  is  proposed  to  reduce  the  cost  evaluations  by  gating 
the  tracks  and  identifying  single  target  and  multiple  target  groupings.  The  tracking  of 
these  groups  is  performed  with  either  a  single  scan  (2D  single-target)  or  multiscan  (ND 
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multitarget)  multidimensional  assignment.  An  overview  of  Algorithm  II  of  MDA  is 
shown  in  Figure  4.  When  multiple  tracks  are  within  a  gate,  the  algorithm  reduces  the 
number  of  cost  evaluations  by  employing  a  filtering  of  the  measurements  by  a  2D  LP 
assignment  via  a  forward-looking  filter  along  the  batch.  These  filtered  measurements  are 
reprocessed  by  an  ND  LP  assignment  MDA.  When  tracks  are  outside  the  vicinity  of  any 
other  tracks — basically,  a  single  target  track — the  algorithm  reverts  to  a  2D  algorithm  as 
outlined  in  Algorithm  I,  with  N  equal  to  2. 

The  2D  assignment  (single-target  mode)  is  very  efficient  due  to  this  linearity  in 
the  number  of  measurements  to  tracks.  For  example,  if  400  scans  were  processed  with 
two  tracks  and  three  measurements  per  scan,  a  total  of  2400  cost  evaluations  or  Kalman 
filters  are  computed  for  this  scenario.  In  the  4D  example  given  in  the  Algorithm  I  section 
of  the  thesis,  with  two  tracks,  three  measurements  per  scan,  with  400  scans  to  be 
processed,  a  total  of  21,600  cost  decisions  or  31,200  Kalman  filters  need  to  be  calculated 
for  this  example.  In  this  case,  the  4D  version  is  13  times  more  expensive  than  the  2D 
version  with  three  measurements  per  scan.  This  savings  is  the  rationale  for  using  the  2D 
LP  when  in  single  target  mode.  When  there  is  no  contention  among  close-by  tracks, 
which  is  the  single  target  case,  a  2D  assignment  is  a  cost-efficient  tracker.  It  is  also  noted 
that  the  2D  assignment  degrades  in  perfonnance  as  the  clutter  level  increases.  In 
simulation,  the  2D  assignment  also  tended  to  exhibit  the  track  repulsion  behavior  for 
long-duration  crossing  targets  more  than  versions  with  longer  batch  lengths. 

It  is  emphasized  that  in  the  multiscan,  multitarget  (ND)  mode  in  Algorithm  II,  the 
Kalman  filter  processes  the  entire  batch  and  then  skips  forward  to  process  the  next  batch 
of  data  out  of  the  gate.  This  further  increases  the  speed  of  the  algorithm  over  the 
processing  from  scan  to  scan  described  for  Algorithm  I.  This  part  of  the  processing  is 
described  in  Section  D.5  for  Algorithm  II. 
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Figure  4. 


Overview  of  Algorithm  II:  2D-ND  LP  MDA  tracking  algorithm. 
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2. 


Track  Initialization 


The  track  initialization  functionality  is  the  same  as  described  in  Algorithm  I.  This 
function  is  performed  at  the  first  scan,  and  states  are  projected  forward  for  the  next  scan. 


3.  Track  Gating 


When  multiple  tracks  are  within  a  gate,  there  are  two  possible  outcomes  from 
Algorithm  II.  When  the  tracks  are  in  gate  proximity,  the  2D-ND  LP  tracker  (multitarget 
tracking)  is  employed  to  process  the  multiple  tracks.  When  the  tracks  are  separated, 
outside  any  tracker  gates,  the  2D  LP  tracker  (single-target  tracking)  is  performed.  The 
gating  procedure  can  be  perfonned  either  spatially  or  temporally.  For  example,  a  spatial 
gating  can  be  performed  between  multiple  tracks  using  the  Kalman  filter  negative  log 
likelihoods  described  in  Algorithm  I  with  an  appropriately  set  gating  threshold.  In  this 
research,  a  time-dependent  gating  procedure  is  used.  Because  the  goal  of  this  study  is  the 
long-duration  crossing,  a  time  window  is  identified  in  which  the  2D-ND  LP  tracker  is 
processed.  The  2D  LP  tracker  is  used  for  all  scans  before  and  after  the  window,  which 
begins  with  time  scan  start  and  ends  with  time  scan  t„  end- 

The  pseudo  code  for  the  2D  LP  tracker  (single-target  case)  is  as  follows: 


if  t„ 


<  K  star t  0r  K 


>t„ 


end  ' 


The  pseudo  code  for  the  2D-ND  LP  tracker  (multitarget  case)  is  as  follows: 
else  (t>t  ,  ,  and  t  <  t  , ). 

\  n  n  start  n  n  end  J 


A  diagram  of  the  time-dependent  track  gating  is  given  in  Figure  5.  This 
methodology  is  used  for  consistency  in  the  simulated  results  given  in  Chapter  III.  A 
spatial  gate  can  be  easily  interchanged  as  part  of  the  gating  procedure  within  this 
paradigm.  Tracks  that  are  about  to  cross  can  also  be  identified  in  time  via  extrapolation 
using  a  linear  least  squares  fit  along  the  track  or  another  procedure.  This  research  uses 
fixed  time  scans  for  switching  between  the  two  paths,  as  shown  in  the  flow  diagram  in 
Figure  4. 
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Figure  5.  Track  gate  switching  is  performed  at  time  scans  tn  start  and  tn  crK|.  The  2D  LP 
tracker  (single-target,  single-scan)  is  used  prior  to  t„_s tart  and  after  t„_end-  The 
2D-ND  LP  tracker  (multitarget,  batch)  is  used  between  time  scans  tn _start  and  tn  end- 


4.  2D  LP  Assignment  Single  Scan 

When  tracks  are  outside  of  a  gate,  a  2D  LP  is  used  to  update  the  state  estimates 
and  the  next  scan  is  processed.  The  costs  are  calculated  as  in  Algorithm  I  outlined  in  the 
Cost  Generation  section.  The  2D  assignment  is  performed  as  outlined  in  the  ND 
assignment  section  of  Algorithm  I  with  N  equal  to  2.  The  modules  GenCost_2D  and 
2DLP_Assign  shown  in  Figure  4  are  special  situations  of  the  generic  GenCost  and 
LP  Assign  modules  as  described  in  Algorithm  I.  Note  that  in  Tables  2  and  5  that  the 
number  of  cost  evaluations  and/or  Kalman  filters  to  process  is  linear  with  the  number  of 
measurements  M  (including  missing  measurement)  in  the  scan  for  the  2D  assignment. 

5.  2D-ND  LP  Assignment  Multiscan 

When  tracks  are  inside  of  a  gate,  a  2D-ND  LP  is  processed.  The  modules  shown 
in  Figure  4  labeled  GenCost_2D_FWD,  2DLP_Assign_FWD,  GenCost  ND,  and 
ND  LP  Assign  account  for  the  main  processing  chain  of  this  algorithm. 
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In  the  GenCost_2D_FWD,  and  2DLP_Assign_FWD,  a  2D  LP  filter  is  run  (A-  1) 
scans  forward  in  time.  This  process  is  similar  to  the  section  in  Algorithm  II  for  the  2D 
LP  Assignment  for  a  single  scan.  The  measurements  that  were  assigned  to  the  track  are 
saved  to  be  reprocessed  by  ND  LP  assignment.  During  the  assignment  process,  it  is 
noted  that  a  track  received  an  actual  measurement  or  a  missing  measurement.  By 
following  this  step,  the  number  of  measurements  mn  processed  by  the  ND  Assignment 

becomes  equal  to  or  less  than  the  number  of  tracks  L  .  In  most  circumstances,  the  number  of 
tracks  is  less  than  the  number  of  measurements.  Note  that  using  the  2D  assignment  as  a  filter, 
there  are  certain  advantages  especially  when  processing  in  low  clutter.  First,  if  all  the  tracks 
are  assigned  an  actual  measurement  on  a  particular  scan  in  the  batch,  this  is  duly  noted  and 
the  missing  measurement  need  not  be  processed  for  that  scan.  Second,  if  none  of  the  tracks  is 
assigned  a  real  measurement,  then  only  the  missing  measurement  needs  to  be  processed  for 
the  scan.  Third,  as  the  clutter  level  A  increases,  the  Kalman  filter  calculation  costs  are 
absorbed  by  the  2D  assignment,  which  is  linear  with  the  number  of  measurements.  The 
filtered  data  sent  to  the  ND  Assignment  are  always  less  than  or  equal  to  the  number  of  tracks 
no  matter  how  large  the  clutter  level.  This  property  has  benefits  and  disadvantages.  In  dense 
clutter,  this  part  of  the  algorithm  is  very  fast  but  performance  also  degrades.  Note  that  the  2D 
assignment  process  discussed  in  this  section  is  only  used  as  a  filtering  process.  The  state 
estimates  derived  from  this  process  are  not  used  at  the  end  of  each  batch.  The  2D  assignment 
filters  the  batch  of  data  and  forms  a  new  filtered  batch  of  data  to  be  reprocessed  by  the  ND 
assignment. 

In  the  GenCost_ND,  and  ND_LP_Assign  sections,  these  modules  are  similar  to 
the  GenCost  and  LP  Assign  shown  in  Figure  1  of  Algorithm  I.  The  difference  is  in 
Algorithm  I,  the  missing  measurement  is  accounted  on  every  scan  of  the  batch.  In  the 
ND  assignment  in  Algorithm  II,  the  2D  assignment  provides  an  index  if  a  missing 
measurement  or  a  real  measurement  is  processed  for  a  particular  scan.  Therefore,  the 
missing  measurement  need  not  be  processed  for  every  scan  in  the  batch.  By  filtering  with 
the  2D  assignment,  the  number  of  2D  LP  filtered  measurements  /„,  to  process  with  the 
ND  assignment  is  less  than  or  equal  to  the  number  of  tracks  L.  The  ND  assignment 
formulation  for  Algorithm  II  has  the  same  set  of  Equations  (25)  and  (26)  as  in  Algorithm  I. 
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The  objective  function  in  Equation  (25)  is  minimized  subject  to  the  constraints  in 
Equation  (26).  The  binary  value  un  is  the  measurement  index  value  filtered  by  the  2D  LP 

forward  fdter.  If  un  equals  0,  a  missing  measurement  is  processed;  otherwise,  the  index  is 
1 ,  representing  the  first  real  measurement.  The  value  mn  equal  to  /„,  is  the  number  of  2D 
LP  filtered  measurements  for  scan  n  from  I  to  (N  -  1 ).  The  measurement  index  i  ranges 
from  0 — the  missing  measurement — to  ( mn  =  /„),  where  l„  is  the  number  of 
measurements  filtered  by  the  2D  forward  filter  per  scan,  and  n  is  the  scan  index. 

The  maximum  number  of  costs  calculated  per  track  for  the  2D-ND  LP  MDA  for  a 
fixed  number  of  measurements  per  scan  in  a  batch  is  given  by  Equation  (27): 

Number  of  costs  (  c  )  =  M(N  - 1)  +  LN~l ,  (27) 

where  M  is  the  number  of  measurements  per  scan  including  the  missing  measurement, 
(A-  1)  is  the  batch  length,  and  L  is  the  number  of  tracks.  The  total  number  of  costs  is 
given  by  Equation  (28): 

Total  number  of  costs  =  (M (N  -1)  +  LN~l)L  .  (28) 

The  total  number  of  Kalman  filters  needed  to  process  a  single  track  in  a  batch  is 
given  by  Equation  (29): 

JV— 1 

Number  of  KF  in  reduced  set  =  M(N  — 1)  +  , 

n= 1 

or  (29) 

LN  -L 

Number  of  KF  in  reduced  set  =  M(N-  1)H - , 

L- 1 

where  L  >  1  in  the  generation  function. 

The  generic  equations  used  to  calculate  the  number  of  costs  and  the  number  of 
Kalman  filters  in  the  reduced  set  for  any  number  of  measurements  in  a  batch  for  a  single 
track  are  given  as  Equations  (30)  and  (31),  respectively. 

JV— 1  JV— 1 

Number  of  costs  (  c  )  =  Z  m„+n  ln’  (30) 

n= 1  n= 1 
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(31) 


JV— 1  N  k- 1 

Number  of  reduced  KF  =  I  mn  +  zip.. 

n= 1  k=2  n—l 

where  ln  is  the  number  of  2D  LP  extracted  measurements  for  scan  n  from  1  to  ( k  —  1 ) 
and  k  is  the  dimensional  assignment  from  2  to  N.  Note  that  ln  <  L  where  L  is  the 
number  of  tracks  for  any  scan  n  . 

6.  Kalman  Filter  Update  Multiscan 

A  Kalman  filter  update  is  perfonned  for  every  scan  in  the  batch  based  on  the 
measurements  obtained  from  the  ND  assignment.  Note  that  this  is  a  different  operation 
than  Algorithm  I,  where  a  sliding  batch  is  employed  and  only  the  first  value  in  the  batch 
is  used  to  update  the  Kalman  filter.  In  Algorithm  II,  the  Kalman  filter  processes  the 
entire  batch,  and  skips  forward  to  process  the  next  batch  of  data  until  the  track  is  out  of 
the  gate.  The  Kalman  filter  follows  Equations  (5)  through  (13)  for  scan  index  ts . 

7.  Example 

As  an  example,  consider  a  case  with  500  scans,  two  tracks  (L  =  2),  and  three 
measurements  including  the  missing  measurement  (M=  3)  that  uses  a  6D  assignment 
(N=  6).  The  tracks  are  separated  during  the  first  200  scans,  crossing  during  the  next  100 
scans,  and  separated  again  during  the  last  200  scans.  A  diagram  of  this  example  is  shown 
in  Figure  6.  Using  Algorithm  I,  the  total  number  of  Kalman  filters  (KF)  for  the  separated 
tracks,  i.e.,  for  400  scans,  in  this  case  is  calculated  as  Num_KF_M_A  *  L  *  Num_scan, 
where  Num_KF_M_A  is  the  number  of  KFs  from  Table  5  for  given  values  of  M  and  N, 
and  Num  scan  is  2  *  200  scans,  or  400.  Thus,  the  total  number  of  KFs  for  the  separated 
tracks  is  363  *  2  *  400,  or  290,400  KF.  From  Table  2,  the  total  number  of  cost 
evaluations  is  calculated  as  cost _M_N,  where  cost _M_N  is  the  number  of  costs  from 
Table  2  for  given  values  of  M  and  N.  Thus,  cost _M_N=  c  *  L  *  Num_scan  =  243  *  2  * 
400,  or  194,400  cost  evaluations. 
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Example  with  Batch  Length  =  5(6D) 


Figure  6.  Example  of  6D  assignment  for  two-target  case  with  500  scans 

and  target  crossings  for  100  scans. 


Using  Algorithm  II  for  the  single-target  section,  where  tracks  are  outside  the  gate 
of  Figure  4  (i.e.,  separated),  the  total  number  of  KFs  for  the  separated  tracks  is  calculated 
as  Num_KF_2D_M_A  *  L  *  Num_scan,  where  Num_KF_2D_M_A  is  the  number  of  KFs 
from  Table  5  for  given  values  of  M  and  N,  and  Num  scan  is  400.  Thus,  the  total  number 
of  KFs  for  the  separated  tracks  using  Algorithm  II  is  3  *  2  *  400,  or  2400  KFs.  From 
Table  2,  the  total  number  of  cost  evaluations  is  calculated  as  cost_2D_M_A,  where 
cost_2D_M_A  is  the  number  of  costs  from  Table  2  for  given  values  of  M  and  N.  Thus, 
cost_2D_M_A  =  c  *  L  *  Num_scan  =  is  3  *  2  *  400,  or  2400  cost  evaluations.  Thus, 
based  on  the  ratio  of  the  total  number  of  KFs,  Algorithm  II  is  121  times  faster  than 
Algorithm  I  for  this  case. 

Using  Algorithm  I  for  the  multitarget  case,  the  tracks  are  within  a  gate  between 
scans  200  and  300,  i.e.,  Num  scan  =  100.  From  Table  5,  the  number  of  Kalman  filters  to 
process  for  Algorithm  I  is  Num_KF_M_A  *  L  *  Num_scan  =  363  *  2  *  100,  or  72,600 
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KF.  From  Table  2,  the  total  number  of  cost  evaluations  cost _M_N  is  calculated  as  c  *  L  * 
Num_scan  =  243  *  2  *  100,  or  48,600  cost  evaluations. 

Using  Algorithm  II  for  the  multitarget  case  requires  two  steps  to  calculate  the 
number  of  Kalman  filters.  The  2D  LP  function  is  used  for  forward  filtering,  and  then  the 
ND  LP  function  is  used  to  calculate  the  cost  based  on  the  2D  LP  filtered  output.  Thus, 
the  total  number  of  KFs  for  the  multitarget  case  =  (Num_KF_2  D_M_N  *  L  *  Num_scan) 
+  (Num_KF  6D_L_/V  *  L  *  Num_scan/(A-  1))  =  (3  *  2  *  100)  +  (62  *  2  *  20),  or  3080 
KFs.  Thus,  based  on  the  ratio  of  the  total  number  of  KFs,  Algorithm  II  is  about  23  times 
faster  than  Algorithm  I  for  the  multitarget  case. 

From  Table  2,  the  total  number  of  cost  evaluations  c  is  calculated  as 
(cost_2D_M_A  *  L  *  Num_scan)  +  (cost_6D_Z_A  *  L  *  Num_scan/(A-  1))  =  (3  *  2  * 
100)  +  (32  *  2  *  20)  =  1880  cost  evaluations.  The  number  of  cost  evaluations  is  reduced 
by  more  than  25  times  by  using  Algorithm  II  instead  of  Algorithm  I.  It  is  important  to 
note  that  much  of  the  savings  is  attributed  to  the  elimination  of  the  sliding  batch  in 
Algorithm  II  and  the  2D  LP  data  filtering  that  reduces  the  number  of  measurements  per 
scan  from  M  to  the  number  of  tracks  L.  The  2D  LP  for  the  single-target  section  of  the 
algorithm  provides  most  of  the  computational  savings.  The  comparison  of  costs  and 
number  of  reduced  Kalman  filters  between  Algorithm  I  and  Algorithm  II  is  provided  in 
Table  7. 


Algorithm  I:  ND  LP 

Algorithm  II:  2D-ND  LP 

N- 1 

N- 1  N- 1 

No.  of  cost  evaluations  c 

Um 

n= 1 

2 x  n; 

n= 1  n= 1 

N  k- 1 

N- 1  N  k- 1 

No.  of  reduced  KF 

ZFK 

k=2  n= 1 

n= 1  k= 2  n= 1 

Table  7.  Comparison  of  Algorithm  I  and  II  functions  for  calculating  the  number  of 
cost  evaluations  and  the  number  of  Kalman  filters  in  a  reduced  set  for  a 
single  track  and  a  batch  length  ( N  -  1),  where  m  is  the  varying  number  of 
measurements  for  scan  n  from  1  to  (k  -  1),  and  k  is  the  dimensional 
assignment  from  2  to  N.  The  value  !„  is  the  number  of  2D  LP  extracted 
measurements  for  scan  n  from  I  to  (k  -  1),  and  k  is  the  dimensional 
assignment  from  2  to  N.  Note  that,  for  any  scan  n,  ln<L  where  L  is  the 

number  of  tracks. 
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III.  SIMULATION  AND  RESULTS 


The  following  simulation  is  based  on  a  study  conducted  by  Willet  et  al.  [1]  for 
various  MHT  batch-style  algorithms  for  a  passive  sonar  line  array.  The  results  from  [1] 
are  not  reproduced  as  part  of  this  thesis.  Figure  7  shows  a  sample  case  with  two  targets 
crossing  in  a  cluttered  environment  where  the  crossing  occurs  over  100  scans.  The  solid 
green  lines  represent  the  true  target  trajectories.  The  zoom  view  on  the  right  shows  the 
difficulty  a  tracker  would  have  in  interpreting  the  data  during  the  long-duration  crossing 
of  100  scans.  The  goal  of  this  study  is  to  extend  the  batch  length,  using  MDA  or  another 
method,  to  beyond  the  length  of  the  crossing.  The  objective  for  the  tracking  algorithm  is 
to  maintain  track  on  all  targets  through  the  crossing.  Because  of  the  exponential  nature  of 
MDA  and  the  current  state  of  computing  resources,  it  is  impractical  to  develop  a  tracker 
with  acceptable  runtime  performance  for  batch  lengths  of  100  scans  or  more  without 
major  pruning  of  the  cost  evaluations  or  the  use  of  other  suboptimal  techniques. 
Algorithm  II  provides  a  more  efficient  method  for  extending  the  batch  length  while  still 
adhering  to  the  basic  concepts  of  MDA. 


time(scans)  time(scans) 


Figure  7.  Sample  simulated  case  with  two  targets  in  linear  straight-line  motion  through 
a  100-scan  crossing  in  a  cluttered  environment,  X  =  0.01  (left).  In  a  zoom  view  in 
the  area  of  the  target  crossing,  or  approximately  100  time  scans  (right),  the  solid 
green  lines  represent  the  true  target  trajectories. 
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In  all  of  the  simulations,  the  Pd  =  70%,  and  or  =  2,  where  Pd  is  the  probability  of 
detection,  and  ar  is  the  measurement  error  standard  deviation.  Thus,  when  or  =  2,  the 
targets  are  within  2 ar  of  one  another  during  the  100  scans  of  crossing.  The  Kalman 
filter’s  initial  process  noise  variance  <72f  equals  0.0625.  In  these  simulations,  the  value 
for  (Tq ,  the  process  noise  standard  deviation  for  the  Kalman  filter,  is  set  to  0.0005,  which 
corresponds  to  a  final  position  standard  deviation  of  T<7q(n(n  +  l)(2«  + 1)/ 6)'  2  =  3.23. 

This  value  is  important  because  it  determines  the  stiffness  of  the  track  and/or  the 
randomness  of  the  simulated  measurements  as  input  to  the  tracker.  X  =  0.0,  0.01,  0.02, 
and  0.05  are  the  clutter  density  levels  for  the  trials.  For  instance,  a  clutter  density  of  0.01 
corresponds  to  a  single  clutter  measurement  in  a  linear  space  of  a  100  points.  The 
negative  Kalman  filter  likelihood  function’s  clutter  density  parameter/^  in  Equation  (14) 
is  set  at  0.01  for  all  cases.  There  were  100  trials  conducted  for  each  value  of  the  clutter 
density  level  X  .  Each  trial  consisted  of  500  scans  of  data  with  two  crossing  targets.  The 
tracks  were  initiated  in  all  scenarios  at  10  and  -10  with  zero  velocity.  Figure  8  shows 
sample  clutter  levels  and  track  crossing  scenarios  used  as  input  and  output  for  the  test 
results.  The  results  from  the  100  trials  actually  included  five  track  crossing  scenarios — 
when  all  tracks  are  tracking  {both  tracking),  both  tracks  bounce  {switched),  the  tracks 
coalesce  {coalesced),  only  a  single  track  is  tracking  {one  only),  and  none  of  the  tracks  are 
tracking  the  target  {neither).  The  coalesced  outcome  seldom  occurred  for  the  MDA 
algorithm — only  for  the  high  clutter  case.  This  reflects  the  constraints  imposed  as  part  of 
the  discrete  optimization  problem.  The  thin  green  and  blue  lines  represent  the  truth  of  the 
target  paths.  The  thicker  solid  green  and  blue  lines  represent  the  tracker  outputs. 
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A  =  0.0  and  both  tracking 


A  =  0.02  and  tracks  switched 
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A  =  0.01  and  only  one  tracking 


A  =  0.05  and  neither  tracking 


Figure  8.  Sample  clutter  levels  and  track  crossing  scenarios  used  as  input  and  output  to 
obtain  test  results:  zero  clutter  scenario  with  both  tracks  on  target  (top  left),  a 
clutter  level  of  0.01  and  only  one  track  on  target  (top  right),  a  clutter  level  of  0.02 
and  both  tracks  bounce  (bottom  left),  and  a  clutter  level  of  0.05  and  neither  track 

is  on  target  (bottom  right). 


Figure  9  shows  the  MDA  results  of  using  Algorithm  I,  a  typical  sliding  batch 
implementation  of  the  MDA-style  tracker,  with  batch  lengths  (N  -  1)  of  1,  2,  and  3 
(corresponding  to  2D,  3D,  and  4D  implementations  of  MDA).  Past  the  4D  mark,  the 
number  of  calculations  to  fonn  the  cost  matrix  is  large  (see  Table  5).  For  4D  or  greater,  it 
is  not  possible  to  obtain  the  simulation  runs  in  a  reasonable  length  of  time.  The  results 
show  that  with  no  clutter  (A,  =  0.0),  the  2D  assignment  tracks  both  targets  in  more  than 
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50%  of  the  cases.  The  exact  value  is  56%,  with  95%  confidence  intervals  between 
46.27%  and  65.73%,  based  on  Bernoulli  trials.  For  batch  lengths  of  2  and  3,  Algorithm  I 
tracks  correctly  more  than  70%  of  the  time  in  no  clutter.  For  a  batch  length  of  2,  72%  of 
the  crossings  are  tracked  correctly,  with  95%  confidence  intervals  between  63.2%  and 
80.8%.  For  a  batch  length  of  3,  the  correct  crossings  increase  to  75%,  with  95% 
confidence  intervals  between  66.51%  and  83.49%.  When  X  =  0.01,  the  correct  tracking 
performance  degrades  to  less  than  30%  of  the  time  for  all  batch  lengths.  Table  8  gives 
the  percentage  of  correct  crossings  (i.e.,  both  targets  tracked  through  the  crossing)  in  the 
Algorithm  I  results,  with  95%  confidence  intervals,  for  batch  lengths  of  1,  2,  and  3  with 
X  =  0.0,  0.01,  0.02,  and  0.05  for  each  batch  length.  Increasing  the  batch  length  from  2  to 
3  does  improve  the  performance  slightly  for  X  =  0.01.  As  the  clutter  level  is  increased, 
the  algorithm  performance  decreases  as  expected. 
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Figure  9.  Results  from  Algorithm  I:  ND  LP  MDA  version  for  straight-line  model, 
<7C,  =  0.0005.  For  each  group  of  bars,  the  three  numbers  at  the  bottom  (1,  2,  3) 
represent  the  batch  length  ( N~  1),  with  100  trials  per  bar.  The  groups,  left  to 
right,  represent  clutter  density  parameter  X  at  0.0,  0.01,  0.02,  0.05. 
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Batch 

Length 

A 

Correct 

Crossings 

(%) 

95% 

Lower 

Bound 

95% 

Upper 

Bound 

1 

0 

56 

46.27 

65.73 

1 

0.01 

1 

0 

2.95 

1 

0.02 

0 

0 

0 

1 

0.05 

2 

0 

4.74 

2 

0 

72 

63.2 

80.8 

2 

0.01 

26 

17.4 

34.6 

2 

0.02 

22 

13.88 

30.12 

2 

0.05 

13 

6.41 

19.59 

3 

0 

75 

66.51 

83.49 

3 

0.01 

27 

18.3 

35.7 

3 

0.02 

21 

13.02 

28.98 

3 

0.05 

20 

12.16 

27.84 

Table  8.  Algorithm  I  correct  track  crossings  (both  targets  are  being  tracked)  in  100 
trials  if  aq  =  0.0005,  with  95%  confidence  intervals.  Batch  lengths  are  1,  2, 
and  3  for  X  =  0.0,  0.01,  0.02,  and  0.05. 


Because  of  the  excessive  time  required  to  process  batch  lengths  of  a  higher  order 
with  Algorithm  I,  Algorithm  II  was  developed  to  speed  up  the  processing  time  by 
combining  the  2D  and  ND  versions  of  the  MDA  algorithm.  The  2D  LP  tracker  is  used 
until  the  tracks  are  about  to  cross.  In  all  cases,  the  crossing  occurs  at  scan  200.  The  2D 
LP  tracker  is  again  used  when  the  tracks  exit  the  crossing  at  scan  305.  As  the  track 
crossing  occurs,  a  2D-ND  LP  tracker  is  used  with  various  batch  lengths.  The  results  from 
these  trials  are  shown  in  Figure  10.  For  the  value  of  oq  =  0.0005,  the  algorithm  performs 
very  well  for  batch  lengths  of  2  through  10,  with  99%  tracking  of  both  targets  in 
moderate  clutter  (/.  =  0.01). 

Table  9  gives  the  percentage  of  correct  crossings  (both  targets  tracked  through  the 
crossing)  in  the  Algorithm  II  results,  with  95%  confidence  intervals,  where  oq  =  0.0005, 
for  batch  lengths  of  1,  2,  3,  5,  and  10  with  X  =  0.0,  0.01,  0.02,  and  0.05  for  each  batch 
length.  The  performance  degrades  with  increasing  clutter  as  expected. 


39 


CD 

O) 

CD 

-t— » 

C 

0 

2 

0 

Q. 


140 


120 


100 


1  2  3  5  10  1  2  3  5  10 

Batch  Length  (N-1 ) 


1  2  3  5  10 


Figure  10.  Results  from  Algorithm  II:  2D-ND  MDA  fast  version,  aq  =  0.0005.  For  each 
group  of  bars,  the  five  numbers  at  the  bottom  (1,  2,  3,  5,  10)  represent  the  batch 
length  (. N  -  1),  with  100  trials  per  bar.  The  groups,  left  to  right,  represent  the 
clutter  density  parameter  X  at  0.0,  0.01,  0.02,  0.05. 
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95% 
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95% 
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Bound 

1 

0 

56 

46.27 

65.73 

1 

0.01 

1 

0 

2.95 

1 

0.02 

0 

0 

0 

1 

0.05 

2 

0 

4.74 
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100 

100 

100 

2 

0.01 

99 

97.05 

100 

2 

0.02 

95 

90.73 

99.27 

2 

0.05 

52 

42.21 

61.79 

3 

0 
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100 

100 
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0.01 

99 

97.05 

100 

3 

0.02 

95 

90.73 

99.27 

3 

0.05 

52 

42.21 

61.79 

5 

0 

100 

100 

100 
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0.01 

99 

97.05 

100 

5 

0.02 

95 

90.73 

99.27 

5 

0.05 

54 

44.23 

63.77 

10 

0 

100 

100 

100 

10 

0.01 

99 

97.05 

100 

10 

0.02 

95 

90.73 

99.27 

10 

0.05 

52 

42.21 

61.79 

Table  9.  Algorithm  II  correct  track  crossings  (both  targets  are  being  tracked)  in  100 
trials  if  aq  =  0.0005,  with  95%  confidence  intervals.  Batch  lengths  are  1,  2, 
3,  5,  and  10  for  X  =  0.0,  0.01,  0.02,  and  0.05. 


Because  the  Algorithm  II  results  were  good  with  oq  set  to  0.0005,  aq  was 
increased  to  0.002  for  the  four  clutter  levels  used  in  the  prior  trials.  These  results  are 
shown  in  Figure  11  for  batch  lengths  of  1,  2,  3,  5,  10,  and  15.  The  results  indicate  that 
the  Algorithm  II  tracking  perfonnance  is  better  than  90%  in  terms  of  correct  crossings  for 
batch  lengths  greater  than  1  with  no  clutter. 

Table  10  gives  the  percentage  of  correct  crossings  (both  targets  tracked  through 
the  crossing)  in  the  Algorithm  II  results,  with  95%  confidence  intervals,  where 
<jq  =  0.002,  for  batch  lengths  of  1,  2,  3,  5,  10,  and  15  with  X  =  0.0,  0.01,  0.02,  and  0.05  for 
each  batch  length.  As  the  clutter  level  increases,  the  tracking  perfonnance  decreases 
proportionately  for  batch  lengths  greater  than  1.  At  each  level  of  clutter  above  X  =  0.0, 
the  tracking  perfonnance  increases  as  the  batch  length  increases.  With  aq  =  0.002,  the 
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trend  in  the  results  for  clutter  levels  X  of  0.01,  0.02,  and  0.05  indicates  that  increasing  the 
batch  length  leads  to  a  higher  probability  of  maintaining  track  for  long-duration  crossings 
of  100  scans. 


1  2  3  51015  1  2  3  51015  1  2  3  51015  1  2  3  51015 


Batch  Length  (N-1 ) 

Figure  11.  Results  from  Algorithm  II:  2D-ND  MDA,  <jq  -  0.002.  For  each  group  of 
bars,  the  six  numbers  at  the  bottom  (1,  2,  3,  5,  10,  15)  represent  the  batch  length 
( N~  1),  with  100  trials  per  bar.  The  groups,  left  to  right,  represent  clutter  density 

parameter  X  at  0.0,  0.01,  0.02,  0.05. 
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69 

59.94 

78.06 
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0.02 

51 

41.2 

60.8 

10 

0.05 

41 

31.36 

50.64 

15 

0 

98 

95.26 

100 

15 

0.01 

72 

63.2 

80.8 

15 

0.02 

56 

46.27 

65.73 

15 

0.05 

44 

34.27 

53.73 

Table  10.  Algorithm  II  correct  track  crossings  (both  targets  are  being  tracked)  in  100 
trials  if  aq  =  0.002,  with  95%  confidence  intervals.  Batch  lengths  are  1,  2, 
3,  5,  10,  and  15  for  X  =  0.0,  0.01,  0.02,  and  0.05. 


The  appendix  presents  the  nonnalized  estimation  error  squared  (NEES)  and  the 
root-mean-square  error  (RMSE)  [25,  pp.  234-235]  for  the  average  when  both  targets  are 
being  tracked  correctly  through  the  long-duration  crossing  for  each  of  the  test  cases  as  a 
function  of  batch  length  and  clutter  density  X.  The  95%  confidence  bounds  for  the  two 
dimensions  (position  and  velocity)  are  plotted  and  labeled  on  all  of  the  NEES  plots.  In 
all  cases,  the  95%  chi  square  distribution  lower  bound  and  upper  bound  are  0.0506  and 
7.3778,  respectively. 

Notice  that  the  tracker  performance  degrades  with  increasing  clutter  density  X  in 
both  the  NEES  plots  and  RMSE  curves  for  both  Algorithms  I  and  II  and  also  in  the  case 
where  aq  is  increased.  This  effect  is  due  to  the  nature  of  the  MDA  and  other  non- 
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Bayesian  association  algorithms,  such  as  nearest  neighbor,  strongest  neighbor,  and  track 
split  approaches  [13,  pp.  119-141],  where  a  Kalman  filter  is  used  to  estimate  the  state 
vector.  A  Kalman  filter’s  estimate  does  not  update  its  covariance  update  to  measure  its 
uncertainty  as  to  whether  the  associated  assignment  originated  from  the  target  or  was 
associated  to  clutter. 

In  the  MDA  batch  approach,  all  measurement-to-track  assignments  are  considered 
along  the  batch,  but  based  on  the  constraints  imposed,  the  maximum  likelihood  set  of 
measurements  is  used  to  update  the  Kalman  filters  associated  to  the  prior  tracks.  In  the 
simulation  being  studied  here,  two  tracks  are  initiated,  and  because  a  standard  Kalman 
filter  is  utilized  for  state  estimation,  any  association  errors  due  to  increasing  clutter  will 
detrimentally  affect  the  tracker  performance.  Longer  batch  lengths  for  Algorithm  I 
should  improve  the  perfonnance  as  clutter  increases  because  the  measurement-to-track 
assignment  is  exhaustive  along  the  batch  and  also  a  longer  time  frame  would  be  available 
to  choose  the  correct  path  from  the  hypothesis  tree.  In  Algorithm  II,  even  though  greater 
batch  lengths  are  more  feasible,  the  approach  is  considered  suboptimal  because  a  2D 
forward  filter  is  used  to  prune  the  measurement  set.  Any  errors  in  the  2D  association 
process  will  influence  the  subsequent  ND  assignment  process  and  affect  the  tracker’s 
results. 

In  Bayesian  association  algorithms,  such  as  PDA  and  MHT,  an  exhaustive  set  of 
probabilities  are  calculated  for  every  measurement-to-track  pair.  In  PDA-style 
algorithms  [13,  pp.  119-141],  the  measurement  origin  uncertainty  is  integrated  as  part  of 
the  tracking  filter,  and  the  covariance  estimate  is  duly  updated  to  properly  measure  the 
consistency  of  a  track.  In  MHT  style  algorithms  [11],  [12],  although  a  Kalman  filter  is 
used  to  estimate  the  state  vector,  the  initialization  of  tracks  is  incorporated  as  part  of  the 
tracking  function,  and  every  possible  hypothesis  is  calculated  for  all  scans  in  time, 
including  the  possibility  that  a  measurement  originated  from  a  new  target.  This  is 
considered  an  optimal  approach. 
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IV.  CONCLUSIONS 


This  study  investigated  two  versions  of  a  batch-oriented,  MDA  tracking  algorithm 
to  examine  target  crossings  on  the  order  of  100  scans  in  duration  as  tracked  by  a  passive 
sonar  line  array.  Many  of  the  computational  load  issues  in  the  generation  and  evaluation 
of  the  costs  were  identified  to  fonn  the  objective  function.  Linear  programming 
relaxation  was  used  to  solve  the  assignment  problem.  A  suboptimal  but  faster  version  of 
the  ND  assignment  was  developed  to  observe  the  effect  of  increased  batch  length  on 
tracking  through  crossings.  Batch  lengths  of  up  to  15  scans,  or  16D  assignment,  were 
developed  and  tested  on  data  with  various  levels  of  clutter  for  the  suboptimal  version. 
The  two  algorithms  were  tested  via  100-trial  Monte  Carlo  simulations. 

Algorithm  I,  the  optimal  version,  was  tested  up  to  the  4D  assignment  case.  The 
results  showed  track  repulsion  behavior  dominating  correct  tracking  through  crossings, 
which  was  achieved  in  less  than  30%  of  the  trials;  the  tracks  were  switched  more  than 
60%  of  the  time  in  moderate  clutter  (k  =  0.01).  With  a  batch  length  of  two  and  X  =  0.01, 
tracks  were  switched  70%  of  the  time,  with  95%  confidence  intervals  between  61.02% 
and  78.98%,  based  on  Bernoulli  trials.  With  a  batch  length  of  three  and  X  =  0.01,  tracks 
were  switched  69%  of  the  time,  with  95%  confidence  intervals  between  59.94%  and 
78.06%.  Table  8  gives  the  percentage  of  correct  track  crossings  for  each  case. 

Algorithm  II — suboptimal  but  faster — exhibited  significantly  better  performance: 
the  percentage  of  correct  track  crossings  exceeded  99%  for  batch  lengths  of  two  or  more 
and  X  =  0.01,  using  the  same  process  noise  standard  deviation  as  for  Algorithm  I.  Table  9 
gives  the  percentages  observed  in  each  case  and  the  confidence  intervals.  When  the 
standard  deviation  was  increased  for  the  Algorithm  II  results,  the  trend  showed  the 
number  of  correct  track  crossings  generally  increasing  as  batch  length  increased.  A 
definitive  conclusion  cannot  be  drawn  on  this  point  because  the  spread  in  the  data,  <10%, 
is  not  statistically  significant.  The  true  effect  of  batch  length  on  correct  crossings  can 
only  be  verified  by  testing  with  batches  much  longer  (i.e.,  (N  -  1)  >  100  scans)  than  the 
15  scans  used  here. 
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Figure  12.  Average  NEES  for  Algorithm  I:  2D  LP,  <jq  =  0.0005  . 


Figure  13.  Average  RMSE  for  Algorithm  I:  2D  LP,  <7q  =  0.0005  . 
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Figure  14.  Average  NEES  for  Algorithm  E  3D  LP,  <7q  =  0.0005  . 


Figure  15.  Average  RMSE  for  Algorithm  I:  3D  LP,  <7q  =  0.0005  . 
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Figure  16.  Average  NEES  for  Algorithm  E  4D  LP,  <7q  =  0.0005  . 


Figure  17.  Average  RMSE  for  Algorithm  I:  4D  LP,  <jq  =  0.0005  . 


55 


Figure  18.  Average  NEES  for  Algorithm  II:  2D-3D  LP,  <7q  =  0.0005  . 


Figure  19.  Average  RMSE  for  Algorithm  II:  2D-3D  LP,  oq  =  0.0005  . 
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Figure  20.  Average  NEES  for  Algorithm  II:  2D-4D  LP,  <7q  =  0.0005  . 


Figure  21.  Average  RMSE  for  Algorithm  II:  2D-4D  LP,  oq  =  0.0005  . 
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Figure  22.  Average  NEES  for  Algorithm  IE  2D-6D  LP,  oq  =  0.0005  . 


Figure  23.  Average  RMSE  for  Algorithm  II:  2D-6D  LP,  oq  =  0.0005  . 
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Figure  24.  Average  NEES  for  Algorithm  II:  2D-1  ID  LP,  <Jq  =  0.0005  . 


Figure  25.  Average  RMSE  for  Algorithm  II:  2D-1  ID  LP,  oq  =  0.0005  . 
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Figure  26.  Average  NEES  for  Algorithm  IE  2D-3D  LP,  <7q  =  0.002  . 


Figure  27.  Average  RMSE  for  Algorithm  II:  2D-3D  LP,  <yq  =  0.002  . 
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Figure  28.  Average  NEES  for  Algorithm  IE  2D-4D  LP,  <7q  =  0.002  . 


Figure  29.  Average  RMSE  for  Algorithm  II:  2D-4D  LP,  <yq  =  0.002  . 
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Figure  30.  Average  NEES  for  Algorithm  IE  2D-6D  LP,  <7q  =  0.002  . 


Figure  3 1 .  Average  RMSE  for  Algorithm  II:  2D-6D  LP,  <yq  =  0.002  . 
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Figure  32.  Average  NEES  for  Algorithm  II:  2D-1  ID  LP,  <Jq  =  0.002  . 


Figure  33.  Average  RMSE  for  Algorithm  II:  2D-1  ID  LP,  <jq  =  0.002  . 
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Figure  34.  Average  NEES  for  Algorithm  II:  2D-16D  LP,  oq  =  0.002  . 


Figure  35.  Average  RMSE  for  Algorithm  II:  2D-16D  LP,  <jq  =  0.002  . 
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